firstname_secondname_thirdname_lastname_02.html


#Import Geopandas for geo-visualization
import geopandas as gpd
#Import Pandas for data analysis
import pandas as pd
#Import Seaborn for data visualization
import seaborn as sns
#Import Contextily for basemaps including OpenStreetMap
import contextily as cx
#Import Matplotlib for data visualization
import matplotlib.pyplot as plt
#Import weights for spatial analysis
from pysal.lib import weights
#Import numpy for array operations
import numpy as np
#Import esda for spatial autocorrelation
import esda


 #I use the geopandas read_file function to read the shapefile.
den_haag_neighborhoods = gpd.read_file("data/neighborhoods.shp")
# I check the data by printing the first 5 rows of the data.
den_haag_neighborhoods.head()


# use plot() to display shapefiles
den_haag_neighborhoods.plot()


#Load data using Pandas read_excel function
csb_data = pd.read_excel("data/kwb-2020.xls")
#Use info() to check the data import
csb_data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 17341 entries, 0 to 17340
Columns: 119 entries, gwb_code_10 to ste_oad
dtypes: float64(33), int64(37), object(49)
memory usage: 15.7+ MB


#By selecting the rows where the column gm_naam is equal to 's-Gravenhage' and remove the duplicates by only choosing the "Buurt" rows
csb_data_den_haag = csb_data[(csb_data["gm_naam"] == "'s-Gravenhage") & (csb_data["recs"] == "Buurt")]
#Afterwards, I will check the first 5 rows of the data to see if the data was loaded correctly using the head() method.
csb_data_den_haag.head()


# display the number of rows and columns
csb_data_den_haag.shape

(114, 119)


# List of columns to keep
measures = ['a_inw','g_woz','a_opl_lg','bev_dich','regio']
#dictionary of columns to rename
measure_dictionary = {
    'a_inw':'Amount of Inhabitants'
    ,'g_woz':'Real Estate Value [x 1,000 €]',
    'a_opl_lg' : 'Amount of people with low education',
    'bev_dich':'Population density [#/km^2]',
    'people with low education [%]' : 'people with low education [%]'
}


#Creating subset with the measure list
csb_data_den_haag_measures = csb_data_den_haag[measures]
#Displaying overview of new subdataset
csb_data_den_haag_measures.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 114 entries, 7766 to 7922
Data columns (total 5 columns):
 #   Column    Non-Null Count  Dtype 
---  ------    --------------  ----- 
 0   a_inw     114 non-null    int64 
 1   g_woz     114 non-null    object
 2   a_opl_lg  114 non-null    object
 3   bev_dich  114 non-null    object
 4   regio     114 non-null    object
dtypes: int64(1), object(4)
memory usage: 5.3+ KB


#Display the first 5 rows of the data
csb_data_den_haag_measures.head()


#Set index to name of the neighborhood
csb_data_den_haag_measures.set_index('regio', inplace=True)
#Display the first 5 rows of the data
csb_data_den_haag_measures.head()


#Display all rows where the column "a_opl_lg" is equal to '.'
csb_data_den_haag_measures[(csb_data_den_haag_measures['a_opl_lg'] == '.')]


#drop neighborhoods with no inhabitants
csb_data_den_haag_measures = csb_data_den_haag_measures[csb_data_den_haag_measures['a_inw']!=0]
#check for remaining rows with missing values
csb_data_den_haag_measures[(csb_data_den_haag_measures['a_opl_lg'] == '.')|(csb_data_den_haag_measures['g_woz'] == '.')|(csb_data_den_haag_measures['bev_dich'] == '.')]


#Create matrix of neighborhood weights with Queen Approach
w_queen = weights.Queen.from_dataframe(den_haag_neighborhoods,idVariable='neighb_cbs')


# Setting amount of people with low education to 0
csb_data_den_haag_measures.at['Kerketuinen en Zichtenburg','a_opl_lg'] = 0
csb_data_den_haag_measures.at['De Reef','a_opl_lg'] = 0


def get_neighborhood_area_average(data,neighborhood,neighbors,measure):
    '''
    This function calculates the average of the measure for the adjacent neighbors of the neighborhood
    Parameters:
        data: the dataset containing the measure
        neighborhood: the name of the neighborhood
        neighbors: the list of neighbors of the neighborhood
        measure: the name of the measure
    Returns:
        The average of the measure for the adjacent neighbors of the neighborhood
    '''
    if measure =='a_opl_lg':
        #Return average Amount of people with low education % of the neighbors   
        return (data.loc[data.index.isin(neighbors),'a_opl_lg'].astype(float)/ data.loc[data.index.isin(neighbors),'a_inw'].astype(int) * 100).mean()
    else:
        # use mean() to calculate the average of the with .loc[] selected subset    
        return data.loc[data.index.isin(neighbors),measure].astype(int).mean()


csb_data_den_haag_measures.at['Westbroekpark','g_woz']  = get_neighborhood_area_average(csb_data_den_haag_measures,'Westbroekpark',list(w_queen['Westbroekpark']),'g_woz')
csb_data_den_haag_measures[(csb_data_den_haag_measures['a_opl_lg'] == '.')|(csb_data_den_haag_measures['g_woz'] == '.')|(csb_data_den_haag_measures['bev_dich'] == '.')]


csb_data_den_haag_measures.at['De Rivieren','g_woz']  = get_neighborhood_area_average(csb_data_den_haag_measures,'De Rivieren',list(w_queen['De Rivieren']),'g_woz')
csb_data_den_haag_measures[(csb_data_den_haag_measures['a_opl_lg'] == '.')|(csb_data_den_haag_measures['g_woz'] == '.')|(csb_data_den_haag_measures['bev_dich'] == '.')]


#To change the data type i use the astype() function on the respective columns
csb_data_den_haag_measures['a_opl_lg'] = csb_data_den_haag_measures['a_opl_lg'].astype(float)
csb_data_den_haag_measures['g_woz'] = csb_data_den_haag_measures['g_woz'].astype(float)
csb_data_den_haag_measures['bev_dich'] = csb_data_den_haag_measures['bev_dich'].astype(float)


# Calculate Percentage of Low Education
csb_data_den_haag_measures['people with low education [%]'] = csb_data_den_haag_measures['a_opl_lg']/ csb_data_den_haag_measures['a_inw'] * 100


# Calculate Percentage of Low Education people with the average of their adjacent neighborhoods for Kerketuinen en Zichtenburg and De Reef
csb_data_den_haag_measures.at['Kerketuinen en Zichtenburg','people with low education [%]']  = get_neighborhood_area_average(csb_data_den_haag_measures,'Kerketuinen en Zichtenburg',list(w_queen['Kerketuinen en Zichtenburg']),'a_opl_lg')
csb_data_den_haag_measures.at['De Reef','people with low education [%]']  = get_neighborhood_area_average(csb_data_den_haag_measures,'De Reef',list(w_queen['De Reef']),'a_opl_lg')


# I do this by using the rename() method on the dataframe. 
# I use the measure_dictionary to rename the columns. 
# I set inplace=True to change the dataframe directly.
csb_data_den_haag_measures.rename(columns=measure_dictionary, inplace=True)
csb_data_den_haag_measures.head()


def graph_distribution(data,measure,axis):  
    '''
    This function plots the distribution of the measure
    Parameters:
        data: the dataset containing the measure
        measure: the name of the measure
        axis: the axis to plot the distribution
    Returns:
        None
    '''
    #deactive grid
    axis.grid(False)    
    #remove the top and right borders
    sns.despine(top=True,right=True)
    #plot the distribution dsitrubted on 20 bins and with kernel density estimation
    sns.histplot(data[measure], kde=True, bins=20,ax=axis)


# Create 3 subplots and define size of figure
figure, axes = plt.subplots(1, 3, figsize=(15, 5))
# Define measures to plot
measures_to_graph = ['Real Estate Value [x 1,000 €]','people with low education [%]','Population density [#/km^2]']
# Loop over the measures and axes and plot them, ravel returns a flattened array for the loop
# zip() function returns an iterator of tuples based on the iterable object
for measure,ax in zip(measures_to_graph,axes.ravel()):
    #plot the distribution using graph_distribution()   
    graph_distribution(csb_data_den_haag_measures,measure,ax)


#use describe() on dataframe column to get descriptive statistics
csb_data_den_haag_measures['Real Estate Value [x 1,000 €]'].describe()

count    111.000000
mean     305.100129
std      156.067976
min      106.000000
25%      181.500000
50%      284.000000
75%      375.000000
max      740.000000
Name: Real Estate Value [x 1,000 €], dtype: float64


#use describe() on dataframe column to get descriptive statistics
csb_data_den_haag_measures['people with low education [%]'].describe()

count    111.000000
mean      21.043396
std       10.914788
min        0.000000
25%       11.637259
50%       20.689655
75%       29.215544
max       48.484848
Name: people with low education [%], dtype: float64


##use describe() on dataframe column to get descriptive statistics
csb_data_den_haag_measures['Population density [#/km^2]'].describe()

count      111.000000
mean      9193.765766
std       5938.341561
min         24.000000
25%       5049.500000
50%       8772.000000
75%      12548.500000
max      22381.000000
Name: Population density [#/km^2], dtype: float64


def graph_scatterplot(data,measure1,measure2,axis):
    """creates a scatterplot using the provided dataframe and x and y measures
    
    Keyword arguments:
        data -- the dataframe to use
        measure1 -- the x measure name
        measure2 -- the y measure name
        axis -- the axis to plot on

    Return: 
        None
    """
    
    #deactive grid
    axis.grid(False)
    
    #remove the top and right borders
    sns.despine(top=True,right=True)
    #graph scatterplot with provided arguments
    sns.scatterplot(x = data[measure1], y=data[measure2],ax=axis)


# Create Subplots with matplot
f,ax = plt.subplots(1,3,figsize=(15,5))
#graph the scatterplots using the graph_scatterplot function and axes of subplots
graph_scatterplot(csb_data_den_haag_measures,'people with low education [%]','Real Estate Value [x 1,000 €]',ax[0]) 
graph_scatterplot(csb_data_den_haag_measures,'Population density [#/km^2]','Real Estate Value [x 1,000 €]',ax[1])
graph_scatterplot(csb_data_den_haag_measures,'Population density [#/km^2]','people with low education [%]',ax[2])


#Use heatmap to show correlation between measures
sns.heatmap(csb_data_den_haag_measures[measures_to_graph].corr(), annot=True)


#replace Kraayenstein with Kraayenstein en Vroondaal, setting inplace = True to change the dataframe directly
den_haag_neighborhoods.replace('Kraayenstein','Kraayenstein en Vroondaal',inplace=True)
#Check if rename was successful
den_haag_neighborhoods[den_haag_neighborhoods['neighb_cbs']=='Kraayenstein en Vroondaal']


#merge the two dataframes using neighb_cbs and regio as keys
den_haag_neighborhoods_with_measures = den_haag_neighborhoods.merge(csb_data_den_haag_measures, left_on='neighb_cbs',right_on='regio')
#Show overview of new merged dataframe
den_haag_neighborhoods_with_measures.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
Int64Index: 111 entries, 0 to 110
Data columns (total 8 columns):
 #   Column                               Non-Null Count  Dtype   
---  ------                               --------------  -----   
 0   neighb_cbs                           111 non-null    object  
 1   neigb_cijf                           111 non-null    object  
 2   geometry                             111 non-null    geometry
 3   Amount of Inhabitants                111 non-null    int64   
 4   Real Estate Value [x 1,000 €]        111 non-null    float64 
 5   Amount of people with low education  111 non-null    float64 
 6   Population density [#/km^2]          111 non-null    float64 
 7   people with low education [%]        111 non-null    float64 
dtypes: float64(4), geometry(1), int64(1), object(2)
memory usage: 7.8+ KB


def create_choroplath_map(data, measure,axis, title):
    """ Creates a choropleth map of the given measure in the given data.
    
        Parameters:
                data (geodataframe_crs): The data to plot
                measure (str): The measure to plot
                axis (axis): The axis to plot on
                title (str): The title of the plot
        Returns:
                None
         """
    #adapt geodataframes to same coordinate system
    den_haag_neighborhoods_with_measures_crs = data.to_crs(epsg=3857)
    den_haag_neighborhoods_crs = den_haag_neighborhoods.to_crs(epsg=3857)
    #plot neighborhoods as base layer showing dropped neighborhoods in grey
    den_haag_neighborhoods_crs.plot(ax=axis, color='grey', edgecolor='black', linewidth=0.5)
    #add layer with choropleth map
    den_haag_neighborhoods_with_measures_crs.plot(column=measure,alpha=0.5,legend=True, legend_kwds={'orientation': "horizontal"},cmap='plasma',ax=axis)
    #disable grid
    axis.grid(False)
    #Set Title to measure name
    axis.set_title(title)
    #disable axis
    axis.set_axis_off()
    #add street view as layer
    cx.add_basemap(ax=ax, url=cx.providers.CartoDB.Positron)


# Create 3 subplots and define size
figure, axes = plt.subplots(1,3, figsize=(30, 10))
# Create a list of measures to plot
measures_to_graph = ['Real Estate Value [x 1,000 €]','people with low education [%]','Population density [#/km^2]']
# Loop over the measures and plot them using create_chlorpleth_map()
for measure,ax in zip(measures_to_graph,axes.ravel()):    
     create_choroplath_map(den_haag_neighborhoods_with_measures,measure,ax,measure)


#figure and subplots for 3 choroplaths
figure, axes = plt.subplots(1,3, figsize=(30, 10))
# Loop over the measures and plot them
for measure,ax in zip(measures_to_graph,axes.ravel()):
    #create a list of the neighborhoods with the 5 highest values  
    highest_values = den_haag_neighborhoods_with_measures.sort_values(by=measure,ascending=False).head(n=5)
    #create choroplath map of the 5 neighborhoods with the highest values
    create_choroplath_map(highest_values,measure,ax, title = measure + ' 5 highest')


#figure and subplots for 3 choroplaths
figure, axes = plt.subplots(1,3, figsize=(30, 10))
# Loop over the measures and plot them
for measure,ax in zip(measures_to_graph,axes.ravel()):
    #create a list of the neighborhoods with the 5 lowest values  
    highest_values = den_haag_neighborhoods_with_measures.sort_values(by=measure).head(n=5)
    #create choroplath map of the 5 neighborhoods with the lowest values
    create_choroplath_map(highest_values,measure,ax, title = measure + ' 5 lowest')


def create_spatial_weights(data, measure):
        """ Creates a spatial weights matrix of the given measure in the given data.
    
        Parameters:
                data (geodataframe_crs): The data
                measure (string): The measure
        Returns:
                Local Indicators of Spatial Association (LISA) object
         """
        #Create a spatial weights matrix
        w = weights.Queen.from_dataframe(data)
        #Transform
        w.transform = 'r'
        #Spatial Lag
        data['w_'+measure] = weights.lag_spatial(w, data[measure])
        #Standardize
        data[measure+'_std'] = (data[measure] - data[measure].mean()) / data[measure].std()
        #standarized weights
        data['w_'+measure+'_std'] = weights.lag_spatial(w, data[measure+'_std'])
        #Create LISA object
        lisa = esda.Moran_Local(data[measure], w)
        #Significance
        data['significant'] = lisa.p_sim < 0.05
        #Quadrants
        data['quadrant'] = lisa.q
        #return LISA object
        return lisa


def create_lisa_cluster_graph(lisa,geodataframe,measure):
    '''
    Creates a graph of the LISA clusters
    Parameters:
                data (geodataframe_crs): The data
                lisa (LISA object): The LISA object
                measure (string): The measure name
    Returns:
        None
    '''

    
    # Setup the figure and axis
    f, ax = plt.subplots(1, figsize=(7,7))
    # Plot insignificant clusters
    
    den_haag_neighborhoods_crs = den_haag_neighborhoods.to_crs(epsg=3857)
    geodataframe_crs = geodataframe.to_crs(epsg=3857)
    # Break observations into significant or not
    geodataframe_crs['significant'] = lisa.p_sim < 0.05
    # Store the quadrant they belong to
    geodataframe_crs['quadrant'] = lisa.q
    #baselayer in white
    den_haag_neighborhoods_crs.plot(ax=ax, color='white', edgecolor='black', linewidth=0.5)
    #plot insignificant clusters in light-grey
    ns = geodataframe_crs.loc[geodataframe_crs['significant']==False, 'geometry']
    ns.plot(ax=ax, color='#dedede')
    # Plot HH clusters in orange
    hh = geodataframe_crs.loc[(geodataframe_crs['quadrant']==1) & (geodataframe_crs['significant']==True), 'geometry']
    hh.plot(ax=ax, color='#eb8134',label='HH')
    # Plot LL clusters in blue
    ll = geodataframe_crs.loc[(geodataframe_crs['quadrant']==3) & (geodataframe_crs['significant']==True), 'geometry']
    ll.plot(ax=ax, color='#6baed6')
    # Plot LH clusters in light-orange
    lh = geodataframe_crs.loc[(geodataframe_crs['quadrant']==2) & (geodataframe_crs['significant']==True), 'geometry']
    lh.plot(ax=ax, color='#b4d5ea')
    # Plot HL clusters in light-blue
    hl = geodataframe_crs.loc[(geodataframe_crs['quadrant']==4) & (geodataframe_crs['significant']==True), 'geometry']
    hl.plot(ax=ax, color='#f2ae7d')

    cx.add_basemap(ax=ax, url=cx.providers.CartoDB.Positron)
    # Style and draw
    ax.set_title('LISA for '+measure)

    #disable axis
    ax.set_axis_off()
    #Show map
    plt.show()


den_haag_neighborhoods_with_measures_for_lisa = den_haag_neighborhoods_with_measures.copy()
# Create LISA object for people with low education [%] using the create_spatial_weights method passing the data and the measure name
lisa_education = create_spatial_weights(den_haag_neighborhoods_with_measures_for_lisa, 'people with low education [%]')
#Displaying the LISA cluster with my created LISA object
create_lisa_cluster_graph(lisa_education, den_haag_neighborhoods_with_measures,'people with low education [%]')


# Create LISA object for people with low education [%] using the create_spatial_weights method passing the data and the measure name
lisa_real_estate_value = create_spatial_weights(den_haag_neighborhoods_with_measures_for_lisa, 'Real Estate Value [x 1,000 €]')
#Displaying the LISA cluster with my created LISA object
create_lisa_cluster_graph(lisa_real_estate_value, den_haag_neighborhoods_with_measures_for_lisa,'Real Estate Value [x 1,000 €]')


# Create LISA object for people with low education [%] using the create_spatial_weights method passing the data and the measure name
lisa_population_density = create_spatial_weights(den_haag_neighborhoods_with_measures_for_lisa, 'Population density [#/km^2]')
#Displaying the LISA cluster with my created LISA object
create_lisa_cluster_graph(lisa_population_density, den_haag_neighborhoods_with_measures_for_lisa,'Population density [#/km^2]')


#show interactive map of Den Haag
den_haag_neighborhoods.explore()

	neighb_cbs	neigb_cijf	geometry
0	Oostduinen	70 Oostduinen	POLYGON ((4.30290 52.12832, 4.30298 52.12827, ...
1	Belgisch Park	71 Belgisch Park	POLYGON ((4.28056 52.11706, 4.28053 52.11706, ...
2	Westbroekpark	73 Westbroekpark	POLYGON ((4.29171 52.10745, 4.29181 52.10739, ...
3	Duttendel	74 Duttendel	POLYGON ((4.32252 52.10768, 4.32252 52.10766, ...
4	Nassaubuurt	48 Nassaubuurt	POLYGON ((4.31943 52.09247, 4.31943 52.09247, ...

	gwb_code_10	gwb_code_8	regio	gm_naam	recs	gwb_code	ind_wbi	a_inw	a_man	a_vrouw	...	g_afs_kv	g_afs_sc	g_3km_sc	a_opp_ha	a_lan_ha	a_wat_ha	pst_mvp	pst_dekp	ste_mvs	ste_oad
7766	BU05180170	5180170	Oostduinen	's-Gravenhage	Buurt	BU05180170	1	0	0	0	...	.	.	.	318	308	10	2597	5	4	683
7768	BU05180271	5180271	Belgisch Park	's-Gravenhage	Buurt	BU05180271	1	8350	4095	4250	...	0,4	0,5	10,9	106	106	0	2587	2	2	2497
7770	BU05180373	5180373	Westbroekpark	's-Gravenhage	Buurt	BU05180373	1	1070	475	590	...	0,7	0,7	18,5	46	41	5	2597	4	2	2402
7771	BU05180374	5180374	Duttendel	's-Gravenhage	Buurt	BU05180374	1	1100	490	610	...	0,4	0,4	14,3	133	130	3	2597	1	3	1135
7773	BU05180448	5180448	Nassaubuurt	's-Gravenhage	Buurt	BU05180448	1	1620	790	830	...	0,2	0,4	22,0	29	28	0	2596	1	1	3307

	a_inw	g_woz	a_opl_lg	bev_dich	regio
7766	0	.	.	.	Oostduinen
7768	8350	398	1120	7880	Belgisch Park
7770	1070	.	60	2591	Westbroekpark
7771	1100	612	90	844	Duttendel
7773	1620	653	80	5737	Nassaubuurt

	a_inw	g_woz	a_opl_lg	bev_dich
regio
Oostduinen	0	.	.	.
Belgisch Park	8350	398	1120	7880
Westbroekpark	1070	.	60	2591
Duttendel	1100	612	90	844
Nassaubuurt	1620	653	80	5737

	Amount of Inhabitants	Real Estate Value [x 1,000 €]	Amount of people with low education	Population density [#/km^2]	people with low education [%]
regio
Belgisch Park	8350	398.0	1120.0	7880.0	13.413174
Westbroekpark	1070	472.4	60.0	2591.0	5.607477
Duttendel	1100	612.0	90.0	844.0	8.181818
Nassaubuurt	1620	653.0	80.0	5737.0	4.938272
Uilennest	2190	380.0	120.0	7439.0	5.479452

EPA-1316 Introduction to Urban Data Science

Assignment 2: Geographic Visualisation¶

`Instructions`¶

1.1 Submission¶

1.2 How do you convert to HTML?¶

1.3 Learning Objectives¶

`Problem`¶

`Data`¶

More Data Sources¶

In case you get more data as shapefiles, and want to play with projections, a nice guide for it is here ¶

`Tasks`¶

`Start your analysis`¶

Analyzing The Hague

1. Import Libraries and Data

2. Problem Statement

3. Dataset

3.1 Loading data

3.2 Data Cleaning

4. Exploratory Data Analysis

4.1 Distribution

4.2 Correlation

4.3 Creating combined dataset

3.4 Choropleths

4. Selection of Neighborhood of Interest

4.1 Local Indicators of Spatial Association (LISA)

4.2 Selection of neighborhoods

	a_inw	g_woz	a_opl_lg	bev_dich
regio
Oostduinen	0	.	.	.
Kerketuinen en Zichtenburg	105	106	.	133
Vliegeniersbuurt	0	.	.	.
De Reef	170	324	.	336
Tedingerbuurt	0	.	.	.

EPA-1316 Introduction to Urban Data Science

Assignment 2: Geographic Visualisation¶

Instructions¶

1.1 Submission¶

1.2 How do you convert to HTML?¶

1.3 Learning Objectives¶

Problem¶

Data¶

More Data Sources¶

In case you get more data as shapefiles, and want to play with projections, a nice guide for it is here¶

Tasks¶

Start your analysis¶

Analyzing The Hague

1. Import Libraries and Data

2. Problem Statement

3. Dataset

3.1 Loading data

3.2 Data Cleaning

4. Exploratory Data Analysis

4.1 Distribution

4.2 Correlation

4.3 Creating combined dataset

3.4 Choropleths

4. Selection of Neighborhood of Interest

4.1 Local Indicators of Spatial Association (LISA)

4.2 Selection of neighborhoods

`Instructions`¶

`Problem`¶

`Data`¶

In case you get more data as shapefiles, and want to play with projections, a nice guide for it is here ¶

`Tasks`¶

`Start your analysis`¶